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THE  USE  OF  ARRAY  ALGEBRA  IN  TERRAIN  MODELING  PROCEDURES 


INTRODUCTION.  Terrain  modeling  is  currently  accomplished  by  a series  of 
polynomial  least-squares  fits  to  observed  elevation  data.  The  polynomials 
are  then  mathematically  joined  to  provide  a second  series  of  polynomials 
that  join  smoothly  with  their  neighbors.  These  final  polynomials  are  then 
used  to  generate  data  for  computer-driven  plotting  devices. 

The  routine  that  produces  the  initial  set  of  polynomials  from  the 
observation  data  utilizes  the  principle  that  the  best  mathematical  fit  of  a 
model  equation  to  arbitrary  points  in  space  must  minimize  the  square  of  the 
vertical  distance  from  the  data  to  the  surface  described  by  the  derived 
function.  This  technique  involves  the  inversion  of  an  m by  m matrix, 
where  m is  equal  to  the  number  of  terms  in  the  model  equation.  Although 
work  at  ETL  currently  uses  a relatively  small  (four  terms)  model  equation, 
new  developments  will  soon  require  more  complexity,  expanding  the  model 
certainly  to  nine  terms  and  probably  more.  Consequently,  the  inversion 
time  for  the  m by  m matrix  will  grow  geometrically  (see  appendix  D) . As 
a means  of  reducing  some  of  these  cumbersome  and  time-consuming  inversions, 
the  technique  of  array  algebra  is  being  investigated.  Array  algebra,  the 
invention  of  Dr.  Urho  A.  Rauhala,  performs  the  same  functions  as  a least- 
squares  fit  in  fewer  computational  steps.  Where  least-squares  requires 
the  inversion  of  one  m by  m matrix,  for  instance,  array  algebra  requires 
the  inversion  of  an  a by  a matrix  and  a b by  b matrix,  where  a and  b 
are  integer  factors  of  m.  As  is  demonstrated  in  appendix  D,  the  computa- 
tional savings  of  array  algebra  increase  dramatically  with  m. 

This  report  analyzes  the  methods  by  which  array  algebra  can  be 
implemented  into  ETL's  terrain  modeling  procedure,  and  determines  its 
feasibility. 

ARRAY  ALGEBRA  IN  TERRAIN  MODELING. 

The  Terrain  Modeling  Procedure.  Prior  to  making  any  determinations 
regarding  the  feasibility  of  array  algebra's  implementation,  it  is 
necessary  to  understand  the  current  terrain  modeling  procedure,  which  can 
be  ordered  into  five  distinct  steps. 

1.  Data  Collection  - This  step  normally  involves  the  analysis  of 
stereo  photography  of  land  defined  by  one  topographic  map  sheet  scaled 
1:50,000.  This  analysis  yields  2.25  million  elevation  observations.  These 
observations,  however,  are  from  stereo  pairs  of  overlapping  photographs  and 
are  organized  accordingly. 

2.  Mosaic  Routine  - This  procedure  uses  a software  routine  to 
transform  step  1 data  into  a set  of  1.02  million  elevations  on  a single 
plane  congruent  to  the  corresponding  1:50,000  map  sheet. 


3.  Preliminary  Least-Squares  (L-S)  Fits  - This  step  establishes 
p number  of  congruent,  overlapping  rectangles  across  the  data  grid  and 
derives  a series  of  fitted  L-S  polynomials,  one  for  each  rectangle.  The 
surfaces  defined  by  these  functions,  however,  do  not  necessarily  join 
smoothly  where  they  overlap. 

4.  Interpolation  - This  step  uses  the  first  set  of  polynomials  and 
a set  of  weighting  functions  to  produce  a second  set  of  polynomials  which 
defines  a set  of  surfaces  that  exhibit  first  order  continuity  (i.e. 
functions  that  agree  in  slope  at  the  points  where  they  meet) . This  set  of 
surfaces  now  defines  a continuous,  smooth-flowing  representation  of  the 
corresponding  terrain. 

5.  Evaluation  of  the  Final  Polynomials  - This  step  involves  the 
evaluation  of  the  final  polynomials  for  each  (x,  y)  at  which  elevation  data 
is  required.  These  elevations  are  then  used  to  produce  contour  maps,  DTM 
data  bases,  3-dimensional  projections,  or  any  of  a variety  of  topographic 
products. 

Array  algebra  can  be  incorporated  into  this  5-step  procedure  only  in 
step  3,  which  currently  uses  a conventional  least-squares  technique,  and  in 
step  5,  which  currently  uses  a standard  algebraic  solution  process.  The 
remainder  of  this  report,  therefore,  will  discuss  only  these  two  steps. 

Deriving  the  L-S  Polynomials.  Each  L-S  polynomial  essentially  converts 
a series  of  data  points  into  a continuous  surface  representation.  The  best 
representation  exists  when  the  polynomials  minimize  the  sum  of  the  squares 
cf  the  vertical  distance  from  the  observations  to  the  defined  surface.  If 
the  fitted  altitude  f(x,y)  is  defined  as, 

f^i’V  = c0  + clxi  + c2yi  + C3Xiyi 
then  the  vertical  distance  d may  be  represented  as 

di  = zi  ~ f<xi’  yi> 

where  z.  is  the  "i"th  observation.  The  entire  set  of  n number  of 
distances  can  then  be  represented  as 

n 

E d. 

. , i 


n 

= E 
i 


= 1 


(Zi  - c0  - 


cjxi  - 


c2yi  ‘ C3Xiyi) 


and  the  sum  of  the  squares  of  these  distances  as 


n 

E (Zi  - cQ  - c1xi  - c2yi  - C3X.y.)2 
i=l 

9 

Minimization  of  E d7  yields,  in  matrix  notation, 
i=l  1 

CL  = (ATA)_1ATZl 

where  Z,  is  an  n by  1 matrix  of  observations,  A is  an  n by  m 
matrix  of  polynomial  terms,  and  is  an  m by  1 matrix  of  coefficients. 

(A  complete  derivation  and  a numerical  example  are  presented  in  appendix  A. 

The  array  algebra  method  of  analysis  follows  a different  approach. 
Assuming  that  the  elevation  data  points  are  on  an  orthonormal  grid  (ordered 
via  an  x/y  grid),  then  we  select  an  e by  f matrix  (ZA)  °f  n 
observations  on  the  grid  and  attempt  an  array  algebra  polynomial  fit.  Then 
in  matrix  notation 

zA  = xcAYT  (I) 

where  X is  an  e by  a matrix  of  x-direction  parameters,  y is  an 
f by  b matrix  of  y-direction  parameters,  and  C.  is  an  a by  b matrix 
of  coefficients.  If  we  use  the  same  m-term  model  equation  that  was  used 
in  the  least-squares  derivation,  then  a and  b must  be  integer  factors 
of  m.  Solving  (1)  for  CA 

CA  = (XTX)-1XTZAY(YTY)-1 

(A  comolete  derivation  and  a numerical  example  are  presented  in  appendix  B) 

It  may  sometimes  be  desirable  to  choose  ab  = m such  that  a ^ b . 

Such  a choice  increases  the  accuracy  of  the  fit  along  the  directional 
parameter  which  corresponds  with  the  greater  factor,  while  denigrating  the 
fit  along  the  other  axis,  thus  producing  a terrain  description  with 
predictable  unidirectional  distortion.  This  may  be  extremely  useful  in 
saving  computational  time  when  analyzing  terrain  that  exhibits  a 
unidirectional  trend  across  a specified  area. 


n 

E 

i=l 
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1 1 
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The  complexity  and  extreme  length  of  Dr.  Rauhala's  array  algebra 
derivation  has  impeded  investigations  into  its  validity,  a fact  that  has 
hindered  its  acceptance  even  in  the  face  of  a high  correlation  of  empirical 
data,  and  there  were  those  who  would  not  accept  the  equivalence  of 


(ATA)-1ArZL  and  (XTX) _1XTZAY (YTY) ~ 1 . 

However,  Mr.  James  R.  Jancaitis,  Automated  Cartography  Branch,  ETL, 
recently  constructed  a detailed  proof^  of  the  equivalence  of  the  array 
algebra  solution  to  that  of  least-squares,  and  also  developed  techniques 
for  weighting  and  constraining  the  array  algebra  fit. 

Impact  Upon  the  Terrain  Modeling  Procedure.  An  observation  of  the 
numerical  examples  in  appendixes  A and  B,  where  identical  best  fits  are 
determined  for  the  same  data  using  both  the  least-squares  and  array 
algebra  techniques,  coupled  with  Mr.  Jancaitis'  proof  of  least-squares 
array  algebra  equivalency,  justifies  the  statement  that  the  implementation 
of  array  algebra  into  ETL's  current  modeling  software  would  have  no  impact 
upon  the  resulting  terrain  model.  The  only  observable  difference  would 
deal  with  the  computational  efficiency  of  the  polynomial  fitting  software. 

As  best  as  can  be  predicted,  therefore,  array  algebra  is  a feasible 
method  for  fitting  a model  equation  to  a set  of  orthonormally  ordered  data. 

AN  ANALYSIS  OF  THE  COMPUTATIONAL  EFFICIENCY  OF  ARRAY  ALGEBRA . 

The  Dual  Nature  of  the  Polynomial  Fitting  Process.  To  measure  the 
relative  efficiency  of  array  algebra,  the  computational  differences 
between  the  least-squares  and  array  algebra  techniques  must  be  isolated. 
This  requires  a quantitative  comparison  of  the  number  of  multiplications 
and  additions  needed  to  compute  the  coefficient  matrices  Cj  or  CA  when 

(ATA)-'AT  Zl 


(XTX)-1XT  ZA  Y(YTY)~l 


CT  = 


'james  R. 
Available 
September 


Jancaitis,  "Theroetical  Analysis  of  Array  Algebra,"  Paper 
at  USA  Engineer  Topographic  Laboratories,  Fort  Belvoir,  VA, 
1976. 
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An  important  feature  of  this  analysis  is  the  repetitive  aspects  of 
these  computations.  Since  each  or  describes  only  one  polynomial, 

then  the  polynomial  must  be  calculated  p number  of  times  to  describe  an 
entire  1:50,000  map  sheet.  The  data  presented  in  Mr.  Jancaitis'  report^ 
leads  to  the  derivation  of  p as 


P 


916  ' 

INTEGER  ROUNDED  UP  

*4(e+l) 


INTEGER  ROUNDED  UP  

*s(f+nj 


where  e & f are  the  respective  number  of  points  on  the  sides  of  each 
rectangle  of  data.  It  should  be  noted  that  p is  subject  to  variations, 
since  the  number  of  observations  on  a 1:50,000  map  sheet  often  vary 
considerably  from  the  example. 

At  first  glance  it  seems  that  the  number  of  multiplications  and 
additions  needed  to  compute  the  coefficient  matrix  must  be  multiplied 
by  p to  predict  the  number  of  computations  per  map  sheet.  However,  this 
is  not  true.  The  computation  of  (A*A)-1at  essentually  involves  only  the 
matrix  A,  which  is  composed  of  m model  equation  terms  evaluated  at  a 
given  (x,y),  listed  over  all  (x,y)  considered  for  each  fit.  The  values  for 
A,  therefore,  are  dependent  only  upon  the  model  equation  and  the  position 
of  the  observations  with  respect  to  the  origin.  If  a new  origin  is 
selected  for  each  computation  of  C|  and  if  it  is  placed  in  the  same 
relative  position  with  respect  to  the  new  set  of  observations,  then  A is 
constant  throughout  the  entire  polynomial  fitting  process.  Thus,  the 
computation  of  (A^A)“^A^  is  a one-time  procedure,  and  only  the 
multiplication  of  (A^A)"^A^  to  Z^  needs  to  be  repeated  p times.  A 
similar  procedure  is  used  with  the  array  algebra  format,  showing  that  the 
computation  of  (X^X)~'x^  and  Y(yTy)~1  are  one-time  calculations  while 
the  multiplication  of  (X^X)~*X’1  to  7. j and  the  multiplication  of 
(XTX)  X^Z^  to  Y(Y^Y)-^  must  be  repeated  p times  per  map  sheet. 

An  analysis  of  multiplications  and  additions  required  to  fit  a model 
equation  of  m terms  to  a series  of  rectangular  data  grids  of  n points 
yields  the  following  figures  (see  appendix  G for  a complete  derivation): 


James  R.  Jancaitis,  "Modeling  and  Contouring  Irregular  Surfaces  Subject 
to  Constraints,  "USA  Engineer  Topographic  Laboratories,  Fort  Belvoir,  VA, 
ETL-CR-74- 1 9 , AD  A0 10406,  January  1975. 


Where 


CT  = (ArA)  *AT  7.j 
Cj  requires: 

2 

l/6m(9mn  + 3n  + 4m  + 3m  - 1)  one-time  multiplications 

mnp  repetitive  multiplications 
2 

l/6m(9mn  - 3n  + 4m  - 6m  - 4)  one-time  additions 

mp(n-l)  repetitive  additions 

Where  CA  = (XTX)_1XT  ?A  Y(YTY)”'  and  ab  = m 


requires: 

1 /6  [a  (9a  ^n+3  ^ n+4a“+3a- 1 ) + b (9b  ^ri+3  ^"n+4b“+b- 1 1] 

one-time  multiplications 

(an  + ab^n)p  repetitive  multiplications 

1/6  [a(9a  ■^n-3'^n+4a^-6a-4)  + b(9b-  ^n-3-^/rH-4b“-6b-4)  ] 

one-time  additions 

(an  - a ^n  + m \T"  -m)p  repetitive  additions 


There  are  some  situations  in  which  some  computational  steps  may  be  saved. 
For  an  explanation  of  these  cases,  refer  to  appendix  C. 

A Comparison  of  the  Computational  Efficiency  of  Least-Squares  to  that 
of  Array  Algebra. 

An  Analysis  According  to  the  Complexity  of  the  Model  Equation.  An 
analysis  was  conducted  (see  appendix  E)  to  determine  the  computational 
savings  of  array  algebra  as  the  model  equation  increased  in  complexity 
from  4 to  256  terms.  To  simulate  undistorted  terrain  modeling,  the  array 
algebra  parameters  selected  were  a = b = \/ro,  where  m was  a perfect 
square,  and  a < b with  a and  b as  near  m as  possible  when  m was 
not  a perfect  square.  The  square  observation  grid  was  selected  over  a 
range  of  n,  compatible  with  the  complexity  of  the  model  equation,  as 


I 


3 

determined  by  Mr.  Jancaitis.  The  results  of  this  analysis  show  that  the 
advantage  of  array  algebra's  smaller  matrix  inversions  is  diminished  by  the 
fact  that  the  inversion  process  is  a one-time  procedure.  The  number  of 
repetitive  multiplications  for  the  conventional  least-squares  procedure, 
however,  exceeds  those  of  the  array  algebra  procedure  by 

mpn  - ap(n  4-  b\/n) 

= mpn  - apn  - apb-y/rT 

= mpn  - apn  - mpY/n 
1 1 

= mpn(  1 -7=^) 

b ym 

Since  b_^2  and  since\/n>3,  then  mpn  (1  - — - >Jp)  is  always  positive, 
and  therefore  array  algebra  requires  fewer  repetitive  multiplications  than 
conventional  least-squares.  Note  also  that  the  difference  increases  with 
b or  Vn-  A similar  result  is  produced  with  the  number  of  repetitive 
additions . 

Appendix  F.  presents  a detailed  listing  of  computer  time  saved, 
considering  only  multiplication  and  addition  as  variable  factors  in  the 
modeling  software.  Use  of  a CDC  6400  computer  is  assumed,  with  no  time- 
sharing problems.  The  data  from  this  listing  show  that  under  the  current 
terrain  modeling  parameters  used  by  Mr.  Jancaitis^,  where 

m = 4 terms 

a = b = 2 

e = f = n 

n = 169  observations 


The  computer  savings  using  array  algebra  will  amount  to  only  41.11 
seconds  per  1:50,000  map  sheet.  More  significant  savings,  however,  can  be 
realized  where  the  model  equation  is  large.  Large  model  equations  allow 
for  better  terrain  resolution  over  a larger  grid,  enabling  larger  values 
for  n to  be  selected  which  in  turn  lowers  the  number  of  polynomials 
required  to  fit  a 1:50,000  map  sheet.  An  attempt  to  do  this  using  the 


3 
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least-squares  technique  falls  because  the  multiplication  and  add  time 
become  unreasonably  large  as  m and  n Increase.  Array  algebra,  however, 
enables  the  procedure  to  be  conducted  without  unduly  large  Increases  in 
computer  time.  Note  figure  1,  which  exhibits  the  relationship  between  the 
length  of  the  model  equation  and  the  necessary  computer  time  for  the 
required  additions  and  multiplications,  assuming  that  n remains  optimal 
for  each  m.  That  is 


n = ( INTEGER  \/42 . 25m) 2 

which  is  a relationship  that  is  derived  from  information  in  Modeling  and 
Contouring  Irregular  Surfaces  Subject  to  Constraints'^.  This  value  for 
n merely  assures  a consistent  terrain  representation  by  matching  the 
size  of  terrain  represented  by  one  polynomial  with  a sufficiently  large 
model  equation. 

An  Analysis  According  to  the  Array  Algebra  Parameters.  The 
multiplication  of 


T — 1 T T - 1 

(X  X)  X1  to  7a  to  Y(Y* Y)  1 


requires  an+mf  multiplications  and  an  + af  + mf  - m additions 
(see  appendix  C) . Since  ab  = m and  ef  = n,  the  number  of  multiplications 
and  additions  increase  as  a is  increased  or  as  f is  increased.  As  a 
result,  the  optimal  values  for  a and  b,  i.e.  those  values  that  minimize 
the  number  of  multiplications  and  additions,  are  a = 1 and  b = m. 

Likewise,  e and  f seem  optimal  at  e = m and  f = 1.  Remember, 
however,  that  a and  b define  the  number  of  terms  in  the  model  equation 
that  carries  x or  y values  and  that  as  a and  b deviate  from 
a - b='yfi’"  the  model  equation  becomes  x-oriented  or  y-oriented,  producing 
a commensurate  lateral  distortion  of  the  fitted  terrain  representation.  A 
similar  situation  develops  as  e and  f deviate  from  e = f=\/n’.  Since 
e and  f define  the  length  of  the  sides  of  the  fitted  rectangular  grids, 
these  grids  deviate  from  squares  to  thin  rectangles.  Since  the  length  of 
either  the  x or  the  y direction  increases,  there  must  be  a corresponding 
increase  in  the  x or  y terms  of  the  model  equation  to  describe  the 
extended  terrain  area.  Thus,  any  variance  of  e must  be  accompanied  by 
a corresponding  variance  in  a,  which  leads  to  the  unidirection  distortion 
described  above. 


Number  of  terms  in  model  polynomial  (M) 


Figure  1:  Time  required  to  calculate 

initial  polynomials 


A quantitative  analysis  of  this  distortion  can  be  performed 
empirically  if  the  present  modeling  software  is  altered  to  incorporate  the 
array  algebra  technique. 

Appendix  F presents  a detailed  listing  of  computer  time  saved, 
considering  only  multiplications  and  additions  as  variable  factors  on  the 
terrain  modeling  software,  with  several  model  equations  studied  for 
differing  values  of  a and  b.  Use  of  a CT)C  6400  computer  is  still 
assumed,  with  no  time  sharing.  The  data  from  the  listing  supports  the 
above  analysis. 

Using  the  Array  Algebra  Format  to  Evaluate  the  Final  Polynomials.  A 
standard  algebraic  method  is  currently  used  to  evaluate  the  final 
polynomials  for  each  (x,y).  This  procedure  requires 

l,018,000q(q  + 2)  adds  per  map  sheet  and 

2 

2,036,000(q  + 2q  - 1)  multiplications  per  map  sheet. 


where  q is  equal  to  the  exponential  order  of  the  model  equation. 

If,  however,  the  array  algebra  x and  y parameter  matrices  are 
defined  as  in  the  polynomial  fitting  process  and  if  the  coefficients  are 
ordered  in  an  x-y  grid,  then 


2 


A 


T 

xcay 


where  X is  an  e'  x a'  matrix 

Y is  an  f'  x b'  matrix 

C.  is  an  a'  x b'  matrix 
A 


The  number  of  multiplications  required  to  produce  7.^  for  all  (x,y)  are 


p'  [(a ' ) 2 (e ' -1)  + b'(e’-l)2] 


and  the  number  of  additions  are 

p'taVe'-lKa'-l)  + (e’-i  )2(b’-l)  ] 
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1 


where,  as  derived  from  Modeling  and  Contouring  Irregular  Surfaces  Subject 
to  Constraints” 


INTEGER  ROUNDED  UP 


%(e'+l)  / V 


INTEGER  ROUNDED  UP 


A detailed  examination  of  these  formulas  is  presented  in  appendix  G. 

The  computer  times  for  the  multiplications  and  additions  necessary  to 
evaluate  all  of  the  polynomials  over  an  entire  map  sheet  were  analyzed, 
and  a comparison  of  conventional  versus  array  algebra  methods  was  conducted. 
The  conventional  system  indicates  a considerable  rate  of  increase  in  com- 
putational computer  time  as  the  length  of  the  model  equation  is  increased, 
a fact  that  has  contributed  to  the  reluctance  to  use  higher  order  polyno- 
mials in  terrain  modeling.  The  array  algebra  format,  however,  enables 
high  order  polynomials  to  be  used  without  significantly  increasing  the 
evaluation  time.  The  multiplication  and  addition  time,  for  instance,  for 
a 4th  order  polynomial  evaluation  over  all  (x,y)  is  56.19  seconds  and 
for  an  18th  order  polynomial  is  168.84  seconds.  As  a comparison,  the 
corresponding  conventional  times  are  26.59  seconds  and  3,235.13  seconds. 

The  following  graph  (figure  2)  indicates  the  results  of  .the  analysis. 

DISCUSSION.  A major  roadblock  in  using  higher  order  polynomials  for 
terrain  modeling  has  been  the  unreasonable  increase  in  required  computer 
time.  Indeed,  the  change  from  a 1st  order  to  a 15th  order  model  equation 
increases  the  computer  time  over  all  five  terrain  modeling  steps  by  almost 
five  hours  per  map  sheet.  A comparative  change  using  array  algebra  results 
in  an  increase  of  only  9 3/4  minutes.  It  should  be  noted,  however,  that 
high  order  polynomials  may  have  drawbacks  that  are  unrelated  to  the 
increase  in  computer  time.  Array  algebra,  by  allowing  these  drawbacks  to 
be  investigated  without  undue  computer  costs,  provides  an  option  that  is 
not  reasonably  available  with  conventional  least-squares. 

When  the  use  of  high  order  polynomials  is  investigated,  the  software 
required  will  be  lengthy  and  exceedingly  intricate.  Since  conventional 
least-squares  cannot  handle  these  polynomials  within  reasonable  computer 
times,  array  algebra  must  be  used.  The  investigation  of  computer  use  of 
high  order  polynomials  and  array  algebra  simultaneously  would  be  excessively 
complex  since  no  workable  basic  structure  exists  for  either  development. 
Plans  at  the  US  Array  Engineer  Topographic  Laboratories,  however,  call  for 
the  current  terrain  modeling  software  (1st  order  polynomials)  to  be  revised 
to  incorporate  an  array  algebra  solution  technique,  followed  by  efforts  to 


^op.  cit. 
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expand  the  model  equation.  Thus,  the  investigation  of  high  order 
polynomials  will  have  array  algebra  available  as  a tool  rather  than  as  a 
problem. 

With  the  advent  of  the  array  algebra  solution  technique  and  especially 
in  light  of  the  recent  work  by  Jancaitis,  the  conventional  least-squares 
software  has  become  at  least  partially  obsolete  for  this  application. 

Since  the  efficient  use  of  computer  time  is  becoming  increasingly  important 
in  virtually  all  aspects  of  automated  cartography,  array  algebra  may  prove 
to  be  a powerful  tool  in  the  search  for  cost-effective  digital  terrain 
modeling  techniques. 

CONCLUSIONS . It  is  concluded  that: 

1.  Array  algebra  produces  the  same  results  as  the  conventional 
least-squares  method. 

2.  Array  algebra  can  be  weighted  and  constrained . ^ 

3.  Array  algebra  performs  a least-squares  type  polynomial  fit  faster 
than  the  conventional  least-squares  method. 

4.  Acceptance  of  either  x-direction  or  y-direction  lateral  distortion 
decreases  array  algebra  computational  time  commensurate  with  the  degree  of 
distortion . 

5.  Array  algebra  evaluates  the  final  polynomial  faster  than  the 
conventional  algebraic  method  in  all  cases  save  that  of  the  smallest 
possible  model  equation  where  array  algebra  is  30  seconds  slower  per  map 
sheet . 


James  R.  Jancaitis,  "Theoretical  Analysis  of  Array  Algebra,"  paper  available 
at  USA  Engineer  Topographic  Laboratories,  Fort  Belvoir,  VA,  September  1976. 
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THE  DERIVATION  OF  C =(ATA)-1A1Z] 

n 

Given  a set  of  n elevation  observations,  Z z±  , suppose  that  one 
wishes  to  fit  the  polynomial  i=l 

f(xi,yi)  = cQ  + Clx.  + c2y.  + c3x.yi 

as  closely  as  possible  to  these  observations.  An  accepted  "best"  fit 
occurs  when  the  sum  of  the  squares  of  the  vertical  distances  from  the 
observations  to  the  described  surface  is  minimized.  The  distance  from  an 
elevation  point  to  the  described  surface  may  be  defined  as 

d = zt  - f^.y^ 

d2  = (z±  - f(xi,yi))2 

l 

The  set  of  n distance-squares  is  therefore 


S = Z &l  = Z (Z4  - c0  - Clx.  - c2yi  - c3xiyi)2 
i=l  i=l 


S will  be  minimized  when  its  derivative  is  equal  to  zero. 


* 6S  = -2E(z.  - cQ  - cjx.  - c2y.  - c3xiyi)  0 

6c0 


SS 

&Ci  = -2E[(z.  - cQ  - c1xi  - c2yi  - 
6S 

6c2  = -2E[(z.  - cQ  - clX.  - ^y.  - 

5S 

6c3  = -2E[(z.  - cQ  - clXi  - c2yi  - 

c0n  + ci  1 xi  + c2  S yi  + c3  1 xiyi 
c0!Xi  + Cj  z*2t  + c2  tx1y1  + c3  £x|yi 

c0  Z yi  + C1  1 Xi  yi  + c2  E yi  + C3 

c0  Z Xiyi  + C1  * Xiyi  + C2  1 X yi  + 


c3xiyi) (xi) ^ 


0 


c x.y.)(y^)]  -0 
3 l l 1 

c3xiyi) (xiyi) ^ = ° 


Ez  i 


X . z . 
1 1 


x.  y±  = s yizi 


2 2 _ 


E x^y.  = E xiyiz1 


n 

*Hereafter  the  symbol  E indicates  E 

i=l 
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The  left  side  of  equation  (A2)  can  now  be  written  as 


which  in  matrix  notation  is 


ATACl 

The  right  side  of  equation  (A2)  can  now  be  written  as 


Equation  (A2)  can  now  oe  rewritten  as 


AtACl  = atzl 

T -1 

Solving  for  , multiply  the  left  side  by  (A  A) 

(ATA)-1(ATA)CL  = (ATA)-1ATzL 
CL  = (ATA)-1  ATZl 


A NUMERICAL  EXAMPLE  OF  THE  LEAST-SQUARES  TECHNIQUE . Given  nine 
observations  on  a 3 by  3 grid. 


z (0 

,0) 

= 100 

z (0 , 

,D 

= 110 

z (0 , 

,2) 

= 112 

z(l. 

,0) 

= 118 

z ( 1 • 

,D 

= 108 

Z ( 1 . 

,2) 

= 120 

112.  ...  * 120.  . . . * 


114 


z (2 ,0)  = 106 


1 


110.  ...  * 108.  ...  * 104 


The  model  equation  is  chosen  as 


Then : 


f(x,y)  = cQ  + cLx  + cxy  + c3xy 


ZL  = 


100 

1 

0 

0 

110 

1 

0 

1 

112 

1 

0 

o 

118 

1 

1 

0 

108 

A = 

1 

1 

1 

120 

1 

1 

2 

106 

1 

2 

0 

104 

1 

2 

1 

114 

1 

2 

2 

T 

A A = 


9 

9 

9 

9 

9 

15 

9 

15 

9 

9 

15 

15 

9 

15 

15 

25 

(aTa)'1  = 1 


36 


25 

-15 

-15 

9 


-15 

15 

9 

- 9 


-15 
9 
15 
- 9 


,.T,rl.T  _ 1 

( A A ) A -^-g- — 


25 

10 

-5 

10 

-15 

-6 

3 

0 

-15 

0 

15 

-6 

9 

0 

-9 

0 

9 

-9 

-9 

9 


-2 


(ATA)  1atz,  = 


L oc 


1_ 

36 


3788 

48 

168 

= 

-36 

105.222 
1.333 
4 . 667 
- 1.000 


c o = 105.222  ci  = 1 1/3  c2  = 4 2/3  c3  = -1 

z = 105.222  + 1 . 333x  + 4.667y  - xy 

Notice  that  this  result  is  identical  with  that  of  the  array  algebra 
numerical  example  in  appendix  B. 


appendix  b. 


ARRAY  ALGEBRA  POLYNOMIAL  FITTING 


THE  DERIVATION  OF  = (XTX)  1XT7-AY(Y  Y) 

When  the  set  of  elevation  observations 


then , 


and,  given  n=ef  observations,  then  all  z . . in  the  grid  are 


Tt  is,  of  course,  possible  to  generalize  this  derivation  by  selecting  a 
general  model  polynomial. 
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Thus,  choosing 


1 1 J_f (XiYj ) _C00+c10xi+c01yj+CllXiy j+C20Xi+C02Vj 


+c  x^v  +c  x v^+c  x^v^+  +c  (a-1)  (b-1) 

+c21xiyj+c12xiyj  C22xiyj  •••  (a-1) (b-1)  i yj 


then , 


(a-1). (b-1) 


Zij 


k,l=0 


c x^y* 
kl  i J 


and  all  z..  in  the  grid  are 


Where  X is  a e by  a matrix  of  x-direction  parameters,  Y is  a f by  b 
matrix  of  y-direction  parameters,  is  an  a by  b matrix  of 

coefficients,  and  ZA  is  a e by  f matrix  of  elevation  observations. 

Solving  for  C 

A 

Z.  = XC.YT 
A A 

ZY  = X C YTY 
A A 

XTZ.Y  = (XTX)C  (YTY) 

A A 

(XTX)-1  XTZaY  = (XTX)~1(XTX)Ca(YTY)  = Ca(YTY) 

T ~1  T T —1  'r  T i 

(XX)  X ZaY(y‘y)  1 = C. (YTY) (YlY)-1  = C 

A A 

T -1  T T i 

CA  = (X  X)  X ZAY(Y  Y)-1 

Note  that  when  X = Y,  a special  computational  case  exists: 
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APPENDIX  C.  THE  NUMBER  OF  MULTIPLICATIONS  AND  ADDITIONS 
REQUIRED  TO  FIT  A POLYNOMIAL  OF  tn  TERMS 
TO  DATA  SETS  OF  n OBSERVATIONS 


REQUIRED  MULTIPLICATIONS  AND  ADDITIONS  USING  THE  METHOD  OF  LEAST-SQUARES. 
Given  that 


cL  = (ata)_1at  zl 


where 


A is  an  nxm  matrix. 


Z is  an  nxl  matrix. 
L 


Then 


T 

L = A A requires 
m 

nEi  multiplications  and 

i=l 


m 

(n-1)  £i  additions 

i=l 


which  simplify  to 

— nm(m+l)  multiplications  and  — m(n-l)vm+l)  additions. 


P = (L)  ' requires 
1 2 

— m( 4m  +3m-l)  multiplications  (see  appendix  D) 


and 


1 2 
1 m(4m  -3m-l) 
o 


additions 


(see  appendix  D) . 


requires 


T 

Q = PA1 

9 

mn  multiplications  and  mn(m-l)  additions. 

R = (A^A)  ^A^  requires 
1 ^ 

_ m(4m“  + 3m  + 9tnn  + 3n  - 1)  multiplications 

6 

and 

1/2 

— m(4m^  - 6m  + 9mn  - 3n  - 4)  additions. 

6 


T -IT 

In  the  polynomial  fitting  process,  it  is  necessary  to  compute  (A  A)  A 
only  once.  The  multiplication  of  (A^A)~^A^  to  Z^,  however,  must  be 
repeated  as  many  times  as  necessary  to  compute  all  of  the  polynomials 
needed  to  describe  the  terrain. 

This  process  requires 

(mn) (#  of  reps.)  multiplications,  and 

m(n-l)(#  of  reps.)  additions. 


REQUIRED  MULTIPLICATIONS  AND  ADDITIONS  USING  ARRAY  ALGEBRA. 

cA  = (XTX)_1XT  zA  Y(YTY)-1 
where 


X 

is  an 

e x a 

matrix 

Y 

is  an 

f x b 

matrix 

ab  = m(terms  in  the  model  equation) 
ef  = n(number  of  observations) 


and 


and 


X X requires 
i.  ae(a  + 1) 

~ a(e  - 1) (a  + 1) 

(X^X)  * requires 

— a(4a1 2  + 3a  - 1) 

6 

1 a (4a2  - 3a  - 1) 

6 

(X^X)  ' X^  requires 

2 

a e 


and 


multiplications 


additions . 


multiplications 


additions  (see  appendix  D) . 


multiplications  . 


ae(a  - 1)  additions. 

T -1  T 

(X  X)  X therefore  requires 

1 2 

_ a(4a  + 3a  + 9ae  + 3e  - 1)  multiplications 
6 


and 


I 


L) 


L 


.1  a(4a  - 6a  + 9ae  - 3e  - 4) 
6 


additions . 


T — 1 

Y (Y  Y)  likewise  requires 


— b(4b*“  + 3b  4-  9bf  + 3f  - 1)  multiplications 

6 


and 


- b (4b2  - 6b  + 9b f - 3f  - 4) 

6 


additions . 


T -IT 

The  repetitive  computations  involve  the  multiplication  of  (X  X)  X 


to  Z and  the  multiplication  of  (X^X)  ^X^Z  to  Y(Y^Y)  *.  These  steps 
require 


p(an  + mf) 


multiplicat ions 


and 


p(an  - af  + mf  - m) 


additions 


where  p is  equal  to  the  number  of  polynomial  fits  per  map  sheet. 
Where  X = Y, 


T - 1 t - 1 

Y (Y  Y)  1 = X(X  X) 


= ^[X(XTX)"l]Tj 


T — 1 T — 1 T 

which,  since  (X  X)  is  symmetric,  is  equal  to  [(XX)  X ] 


T - 1 

therefore,  Y(Y  Y)  need  not  be  calculated. 
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APPENDIX  D..  THE  NUMBER  OF  MULTIPLICATIONS  AND  ADDITIONS 
NECESSARY  TO  INVERT  AN  mxm  MATRIX 
USING  GAUSSIAN  ELIMINATION 


Given  that: 


10  0 
0 10 
0 0 1 

0 0 0 


X 0 0 
X 1 0 
X 0 1 

X 0 0 


0 

0 

0 

1 

0 

0 

0 


1 

0 

bl3 

• • • blm 

X 

x n ... 

0 

0 

1 

C23 

* # ’ C3m 

X 

X o ... 

0 

0 

0 

C33 

C 3m 

X 

X 1 ... 

0 

C = 

0 

0 

Cm3 

. . . c 

mm 

X 

X o ... 

1 

Gauss i an 

reduction 

f rom 

A to 

B requires 

m“ 

multiplications 

and 

m(m-l) 

additions 

. Reduction 

from  B to 

C 

requires  an  identical 

number  of  operations.  Thus,  a complete  mxm  matrix  inversion  requires 


m 


3 


multiplications 


and 


m^"  (m-1) 


additions . 


When  the  original  matrix  is  symmetric,  however,  procedures  are  available 
that  reduce  the  required  multiplications  and  additions  by  a considerable 
amount . 

Diagonalizing  the  matrix,  we  see  that  reduction  from  A to.  B requires, 
as  above, 


m 


multiplications 


and 


m(m-l) 


additions. 
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The  next  step  in  the  diagonalization  involves  reducing  B to  C, 
where 


C = 


1 hl2  b13  • • • \m 


0 1 c ...  c 

23  3m 


0 0c 
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' 3m 


0 0 cm3  ' ' ' cmm 


X 0 0 


X X 0 


XXI 


X X 0 ...  1 


This  reduction  requires 


m(m-l) 


multiplications 


and 


m(m-2) 


additions . 


The  patterns  now  become  apparent,  with  diagonalization  requiring 


i(m-0)+m(m-l)+m(m-2)+  . . . +m[m-(m-l)l  multiplications 


and 


mj^i-l)+m(m-2)+m(m-3)+  . . . +m[m-(m)]  additions 
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-iW 


or, 


in  summation  notation 


■ 


m-1 

£ m(m-i) 
i-0 


and 


m 

v m(m-i) 

L=1 


turn  simplify  to 


'2m2  (nH-1) 


and 


multiplications 


addi t ions 


multiplications 


2 

’2 m (m-1)  additions. 


will  now  consider  the  Gaussian  reduction  of  the  upper  right  half  of  a 
vmnk  trie  matrix  that  has  been  reduced  to  the  diagonal  matrix  A'  : 


1 

a12  * 

• ai(m-l) 

a 

lm 

w 

11 

0 . . 

. 0 

0 

0 

1 

' a2 (m-1) 

a2m 

w 

21 

W22  • ' 

. 0 

0 

0 

0 

. 1 

a3m 

w(m-l)  1 

W(m-1) 2 • • 

W(m-l)(m-l) 

0 

0 

0 

. 0 

1 

wml 

wm2  ' • 

wm(m-l) 

w 

mm 
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Since  the  upper  right  section  of  the  inverse  (containing  zeros  in  A'  ) 
will  eventually  be  defined  by  the  lower  left  section,  it  is  not  necessary 
to  perform  any  operations  on  these  zeros. 

Thus,  where 


The  operations  converting  A'  to  B'  require 


(m-1)  + (m-2)  + (m-3)  + . . . + [m-(m-l)] 


m-1 
= 2 i 

i=l 


Continuing  this  process  presents  the  series 


m-1  m-2  m-3  m-(m-l) 

2 i + Z i + I i + . . . + Z i 
i=l  i=l  1=1  i=l 


or 


'2  (m-l)m+ '?(m-2)  (m-l)+ i2(m-3,i  (m-2)+  . . . + '2  (m-(m-l)  ] fm-(m-2)  ] 


or 


m-1 

*5  £ (m-i)  [m-(i-l)  ] 

i=l 

which  represents  the  number  of  multiplications  or  additions  required  to 
find  the  upper  diagonal  and  can  be  simplified  to 


■lm(m^-l) 

6 

Total  multiplications  required  to  invert  a symmetric  mxm  matrix  are 
therefore,  equal  to 


%m2  (m+l)dim(m 

6 


-i) 


which  simplifies  to 


Im(4m2+3u-l)  multiplications. 

6 

Total  additions,  likewise,  are  equal  to 


\ m^ (m-l)+im(m2-l) 

6 

which  simplifies  to 


2 

■rm(4m  -3m-l) 
6 


38 


APPENDIX  E.  LEAST-SQUARES  VS.  ARRAY  ALGEBEA:  A 

COMPARISON  OF  COMPUTATIONAL  TIMES 
AS  THEY  VARY  WITH  CHANGES  IN  THE 
LENGTH  OF  THE  MODEL  EQUATION 


A comparison  of  multiplication  and  addition  times  in  programs  of  polynomial 
fitting  using  both  the  least-squares  and  array  algebra  techniques. 

Assuming  the  use  of  a CDC  6400  computer.  Multiplication  time  for  the 
CDC  6400  is  5700  nanoseconds;  addition  time  is  1100  nanoseconds. 
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APPENDIX  F.  LEAST-SQUARES  VS.  ARRAY  ALGEBRA:  A 

COMPARISON  OF  COMPUTATIONAL  TIMES  AS 
THEY  VARY  WITH  THE  ARRAY  ALGEBRA 
DIRECTIONAL  PARAMETERS  (e  AND  f ) 


A comparison  of  multiplication  and  addition  times  in  programs  of  polynomial 
fitting  using  both  the  least-squares  and  array  algebra  techniques. 

Assuming  the  use  of  a CDC  6400  computer.  Multiplication  time  for  the 
CDC  6400  is  5700  nanoseconds;  addition  time  is  1100  nanoseconds. 
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APPENDIX  G: 


EVALUATING  THE  FINAL  POLYNOMIALS. 


The  Conventional  Evaluation.  The  evaluation  of  the  qth  order 
polynomial 


? 2 2 2 

z = Cg+c ^x+C2y+c^x  y+c^xy  +c^x  y + 


+c_xTy^ 


has  previously  been  accomplished  by  first  creating  two  arrays,  one  of 
exponential  values  of  x,  another  of  y.  Creating  these  arrays  required 
2 (q— 1 ) multiplications.  The  actual  evaluation  then  required 


q 

£ ( 2 i + 1)  additions  and 

i=l 


q 

£ ( 4 i ) multiplications. 

i=  1 


Total  adds  and  multiplications  were  therefore 


q(q  + 2)  additions,  and 

2 

2(q  + 2q  - 1)  multiplications. 

This  procedure  was  repeated  for  every  desired  elevation,  normally  1.018 
million  times. 

The  Array  Algebra  Evaluation.  The  array  algebra  format,  however, 
lends  itself  to  a more  efficient  method  of  model  equation  evaluation. 
Consider  the  grid  of  polynomials  in  figure  Gl.  Given  these  polynomials, 
we  are  required  to  produce  elevation  data  at  each  "X".  (Although  figure 
Gl  shows  each  polynomial  describing  n'  = 16  elevation  points,  n'  is 
not  restricted  to  this  value,  nor  is  this  4 by  4 arrangement  necessarily 
optimal).  Now  we  define  the  local  coordinate  system  to  describe  the 
relative  positions  of  the  elevation  data  within  the  polynomial  p.  only, 
(see  figure  G2) . 
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Figure  Gl.  The  Grid  of  Polynomials 


. 14 X 

I (3,0)  (3,1)  (3,2) 

X X 

(2,0)  (2,1)  (2,2) 

X X 

(1,0)  (1,1)  (1,2) 

X X 

(0,0)  (0,1)  (0,2) 


(3.3) 

(2.3) 

(1.3) 
(0,3) 


Figure  G2.  The  Local  Coordinate  System 


If  the  matrix 
model  equation  and 
values  of  the  model 


X is  defined  as  an  array  of  x-direction 
the  matrix  Y is  defined  as  an  array  of 
equation,  and  if  the  model  equation  is 


values  of  the 
y-direction 
defined  as 


c00+c10xi+c20x^cllxiyi+c21x^i+C02y^12Vi+C22Xiyi 


then 


1 0 0 

1 0 0 

1 1 1 

1 1 1 

1 2 4 

Y = 

1 2 4 

1 3 9 

1 3 9 



— — 

The  matrix  of  coefficients  defining  each  polynomial  is,  of  course. 


00 

coi 

n 

O 

ro 

10 

C11 

C 1 2 

:20 

C21 

c?2 

_ 

— 

Note  that  the  number  of  terms  in  the  model  equation  (m')  is  chosen  as  9, 
and  is  an  a'  by  b'  matrix,  where  a'b'  = m.  The  model  equation  can 

now  be  written  in  matrix  notation  as 


ZA  = XCAyT 

where  Z.  is  an  e'  by  f'  matrix  of  elevation  points,  and  e'f'  = n'. 

A 

A re-inspection  of  figure  G1  shows  that  since  many  of  the  Z values 
lie  on  borders  of  two  or  more  polynomials,  not  all  Z values  need  be 
calculated  for  every  polynomial. 
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Figure  G3.  Evaluation  of  the  Grid  of 
Polvnimials 

L Rather  than  e'  by  f'  elevation  points,  only  an  (e'-l)  by  (f'-l)  array 

(see  figure  G3)  need  be  calculated.  Also,  an  inspection  of  the  X and 
Y matrices  shows  that  they  are  independent  of  the  coefficients  of  the 
model  equation  and  dependent  only  upon  the  size  of  the  x-y  positional 
grid.  Since  each  polynomial  has  an  origin  and  a directional  grid  (figure 
G2)  that  is  congruent  to  each  origin  and  grid  in  the  remaining  polynomials, 
then  X and  Y (and  of  course  Y^)  are  constant  for  the  entire  map  sheet 
and  need  be  calculated  only  once.  In  a general  case,  the  number  of 
repetitive  multiplications  for  = XCY^  are 

p'  [(a’)2(e'-l)  + b ' (e'-l)2] 


and  repetitive  additions  are 
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g 

Mr.  Jancaitls  describes  Che  number  of  polynomials  ( p ' ) as  a function 
of  the  number  of  elevations  on  the  sides  of  the  grids  upon  which  the 
final  polynomials  are  defined: 

') 
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